home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / rwvector.lha / RWVector2.1 / src / mathpack / dcffti1.f < prev    next >
Text File  |  1989-08-14  |  2KB  |  76 lines
  1.       subroutine dcffti1 (n,wa,ifac)
  2. c
  3. c     Double precision version.  -tk
  4. c
  5. C***BEGIN PROLOGUE  DCFFTI1
  6. C***REFER TO DCFFTI
  7. C***ROUTINES CALLED  (NONE)
  8. C***END PROLOGUE  DCFFTI1
  9.       implicit double precision (a-h,o-z)
  10.       dimension       wa(1)      ,ifac(1)    ,ntryh(4)
  11.       data ntryh(1),ntryh(2),ntryh(3),ntryh(4)/3,4,2,5/
  12. C***FIRST EXECUTABLE STATEMENT  DCFFTI1
  13.       nl = n
  14.       nf = 0
  15.       j = 0
  16.   101 j = j+1
  17.       if (j-4) 102,102,103
  18.   102 ntry = ntryh(j)
  19.       go to 104
  20.   103 ntry = ntry+2
  21.   104 nq = nl/ntry
  22.       nr = nl-ntry*nq
  23.       if (nr) 101,105,101
  24.   105 nf = nf+1
  25.       ifac(nf+2) = ntry
  26.       nl = nq
  27.       if (ntry .ne. 2) go to 107
  28.       if (nf .eq. 1) go to 107
  29.       do 106 i=2,nf
  30.          ib = nf-i+2
  31.          ifac(ib+2) = ifac(ib+1)
  32.   106 continue
  33.       ifac(3) = 2
  34.   107 if (nl .ne. 1) go to 104
  35.       ifac(1) = n
  36.       ifac(2) = nf
  37.       tpi = 6.28318530717959
  38.       argh = tpi/float(n)
  39.       i = 2
  40.       l1 = 1
  41.       do 110 k1=1,nf
  42.          ip = ifac(k1+2)
  43.          ld = 0
  44.          l2 = l1*ip
  45.          ido = n/l2
  46.          idot = ido+ido+2
  47.          ipm = ip-1
  48.          wldc = 1.
  49.          wlds = 0.
  50.          arg = float(l1)*argh
  51.          dl1c = cos(arg)
  52.          dl1s = sin(arg)
  53.          do 109 j=1,ipm
  54.             i1 = i
  55.             wa(i-1) = 1.
  56.             wa(i) = 0.
  57.             ld = ld+l1
  58.             wldch = wldc
  59.             wldc = dl1c*wldc-dl1s*wlds
  60.             wlds = dl1s*wldch+dl1c*wlds
  61.             dc = wldc
  62.             ds = wlds
  63.             do 108 ii=4,idot,2
  64.                i = i+2
  65.                wa(i-1) = dc*wa(i-3)-ds*wa(i-2)
  66.                wa(i) = ds*wa(i-3)+dc*wa(i-2)
  67.   108       continue
  68.             if (ip .le. 5) go to 109
  69.             wa(i1-1) = wa(i-1)
  70.             wa(i1) = wa(i)
  71.   109    continue
  72.          l1 = l2
  73.   110 continue
  74.       return
  75.       end
  76.